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Abstract 

The effect of columnar pins on the flux-hnes melting transition in high-temperature supercon- 
ductors is studied using Path Integral Monte Carlo simulations. We highlight the similarities and 
differences in the effects of columnar disorder on the melting transition in YBa2Cu307_5 (YBCO) 
and the highly anisotropic Bi2Sr2CaCu208+5 (BSCCO) at magnetic fields such that the mean sep- 
aration between flux-lines is smaller than the penetration length. For pure systems, a flrst order 
transition from a flux-line solid to a liquid phase is seen as the temperature is increased. When 
adding columnar defects to the system, the transition temperature is not affected in both materials 
as long as the strength of an individual columnar defect (expressed as a flux-line defect interaction) 
is less than a certain threshold for a given density of randomly distributed columnar pins. This 
threshold strength is lower for YBCO than for BSCCO. For higher strengths the transition line is 
shifted for both materials towards higher temperatures, and the sharp jump in energy, characteris- 
tic of a flrst order transition, gives way to a smoother and gradual rise of the energy, characteristic 
of a second order transition. Also, when columnar defects are present, the vortex solid phase is 
replaced by a pinned Bose glass phase and this is manifested by a marked decrease in translational 
order and orientational order as measured by the appropriate structure factors. For BSCCO, we 
report an unusual rise of the translational order and the hexatic order just before the melting 
transition. No such rise is observed in YBCO. 
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I. INTRODUCTION 



Type II superconductor3i>2j^ allow for a partial penetration of magnetic field into the bulk 
of the superconducting (SC) material when the applied field H satisfies Hd < H < Hc2- 
In a seminal work Abrikoso'v^ showed that when the ratio A/^, where A is the magnetic 
field penetration depth and ^ is the coherence length, is greater than 1 / the magnetic 
field penetrates the SC material in the form of flux-lines (FLs). These FLs are also called 
vortices, since they are surrounded by circular currents. Each FL carries a quantized unit 
of flux 00 = hc/2e called the fluxoid. The FLs have cylindrical shape of radius ~ A (the 
radius is not sharp since the magnetic field decays exponentially like exp(— r/A), where r is 
the distance from the axis) and a non-SC core of radius ~ C,- Due to a repulsive interaction 
among the Fls, they arrange themselves in a triangular lattice referred to as the vortex solid 
(VS). This result follows from mean- field theory. 

After high-temperature superconductors were discovered in the 1980's, it became appar- 
ent that thermal fluctuations, not included in the mean-field theory^, play an important 
role at relatively high temperatures and fields, still below Tc and ifc2- These fluctuations 
can cause the Abrikosov lattice to melt into a disordered liquid via a first order transition 
(FOT), which can be roughly estimated using the Lindemann criterion known from solid 
state physics^'^'*^. Technologically, the melting of the FL lattice is important since the vor- 
tex liquid (VL) is not actually SC due to the dissipation caused by the FL motion when an 
electric current passes through the system. Pinning of FLs by naturally occurring defects in 
the form of vacancies, interstitials, twin and grain boundaries etc., is effective to impede FL 
motion in the VS phase, where the FLs form a rigid correlated network. The effectiveness 
of the pinning manifests itself by leading to high critical currents. In the VL phase pinning 
of a few vortices does not inhibit others from moving when a current is applied. Thus for 
practical purposes the sudden increase in resistivity occurs at the melting transition rather 
than when H = Hc2{T) for any reasonably non- vanishing currents. 

The existence of the melting transition in high-Tc pristine materials has been established 
through numerous experimental^ i^^i^^i^^'^? and numerica l^^'-^^'^i^^i^'^i^^i^'^i^^ studies. As was men- 
tioned above, disorder in the form of points defects and sometimes more extended defects 
can and does occur naturally in laboratory samples. In addition artificial point defects can 
be induced by bombarding the sample with electrons originating from particle accelerators. 
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Extended columnar defects in the form of linear damaged tracks piercing through the sample 
can be induced by heavy ion irradiation. Both naturally occurring and artificially induced 
defects are situated at random positions in the sample and their effective pinning strength 
(i.e. their interaction with FLs) can also vary from defect to defect. Thus defects play 
the role of quenched disorder. The adjective "quenched" refers to the immobility of these 
defects during experimental time scales. Introduction of disorder in terms of point defects 
or columnar pins affects both the properties of the solid and liquid phases and might also 
shift the location of the melting transition in the H-T plane^*^. In the case of point pins, the 
VS phase is replaced with a Bragg glass phase^i*^, characterized by quasi-long-range order. 
The melting transition is predicted to shift towards lower temperature a^^i^^i^^ . In the case 
of columnar pins the VS phase is replaced with a so called pinned Bose glass^^ where FLs 
are trapped by the columnar defects and the whole lattice becomes immobile. The Bose 
glass phase is similar to the localized phase of a two dimensional repulsive Bose gas in the 
presence of quenched disorder, as will be explained in more detail in the next section. 

The effect of both kinds of disorders on the FLs melting has been studied experimen- 
tally in various high-temperature superconductors. Two common materials that have been 
extensively investigated are YBa2Cu307_5 (YBCO) and Bi2Sr2CaCu208+5 (BSCCO), both 
having critical temperatures ranging between 90-120 K at H = 0. The main difference be- 
tween these materials is their anisotropy parameter 7^ = mz/m± > 1, where and m± 
denote the effective masses of electrons moving along the c-axis and the afe-plane respectively. 
BSCCO is much more anisotropic: its anisotropy parameter 7 lies in the range of 50-200 
compared to the range of 5-7 for YBCO^. This fact causes the FLs to be much "softer" 
or elastic. Thus in the case of BSCCO the FLs are sometimes described as a collection of 
loosely connected "pancakes" residing in adjacent Cu-0 planes. Experimental studies on 
YBCO have shown a marked shift in the irreversibility line in the presence of the columnar 
disordeii^IiSSi^SiSli^. The irreversibility line in the H-T plane marks the onset of hysteresis 
effects and is located close to the melting transition on the solid phase side. For BSCCO, 
many experimental studies have been conducte d^^'^^'^'^'^^i^'^^'^'^^i^'^'^'^'^^i^'^'^ . The more recent ones 
have shown!2&i2L2Si^ that the melting line is not shifted when the density of columnar defects 
is relatively low, « but for B^ B a, shift in the position of the melting transition 
is observed. Here the matching field B^ is defined as B^ = ndcj)^ where is the density of 
the columnar defects and 0o is the flux quantum. 
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Theoretical work on columnar disorder includes Bose Glass theory2&. RadzihovskjM^ 
considered the possibility of two kinds of Bose glass phases (strongly or weakly pinned) 
depending on whether B < ot B > B^. More recently, columnar as well as point disorder 
were investigated by Goldschmidt^ using replica field theory. He showed that the melting 
line shifts to lower temperature in the case of point disorder and to higher temperature in 
the case of columnar disorder. 

Due to the complexity of the problem, especially in the presence of disorder, simulations 
have been very useful in studying the FL system. There have been many simulation studies 
of the vortex system in the presence of disorder. However, most simulation work has been 
confined to the addition of point disorder only. In particular, there is little work done on 
the effects of columnar disorder on the FL melting. Recent work by Wengel and Tauber— 
concentrated on the case of high defect density region B^ ^ B . In contrast, a recent 
simulation study by Sen et ah^"^ uses a small density of columnar defects, each of infinite 
strength, but considers an extremely small magnetic field. Similarly, Nandgaonkar et al^ 
also investigate the case of a very small magnetic field. In this paper we consider columnar 
defects of a finite strength {rf) with a relatively low defect density, B^/B = 0.2, and with 
realistic magnetic fields as used in the experiments. 

II. THE MODEL 

We first discuss the method implemented for YBCO: 

Following Nelson^, we map the system of vortices (FLs) in a high-T^ superconductor 
onto interacting bosons in 2-space + 1-time dimensions. Then we do a Path Integral 
Monte Carlo (PIMC) simulation on this system. The partition function of vortices can 
be expressed as: 



where Ri(-2) denotes the 2D position vector of the i'th vortex at a height z along the c-axis, 
F is the free energy as a function of temperature T, is the length of the sample along 
the z-direction. The London free-energy functional 3 is given by 




(2.1) 




(2.2) 
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where Eq = 0o/ (4vrA)^ is the vortex hne energy per unit length, the hne tension is ei 
£o ln(ao/2y^,^)/7^ and Rij{z) = \Ili{z) — Rj(z)|. Here 



.0^,/^ (2.3, 

is the lattice spacing, B is the magnetic field along the z-direction and 7 is the anisotropy 
parameter. Ko{Rij/X) is the in-plane interaction potential between two FLs a distance 
Rij apart. Kq denotes the modified Bessel function of the first kind. This expression for 
the London free-energy is an approximation that neglects the non-local interaction of real 
vortices and replaces it by in-plane interactions only, which is really justified if the FLs do 
not deviate too much from straight lines along the z-axis^^. With this approximation, the 
system of interacting FLs is equivalent to a system of bosons in d = 2 dimension 
interacting with a pairwise potential Ko{Rij/X). 

The path integral representation of a system of A^ bosons of mass m each in two dimen- 
sions, interacting through a potential V{r) = g'^KQ{r/X), with g being the strength and A 
being the range of the repulsive interaction, is given at finite temperature Tg in terms of 
the imaginary time action 




h h In I ^ 2 \ dr J ^ V A 



(2.4) 



Here r is the imaginary time and Tg is the temperature of the Bose system. We see that 
there is a one to one parameter mapping between the boson system and the vortex systenJ^: 

r ^ z, kT, ^ 2eo, h/kTe ^ L„ ei, n = N/A = 2/{V3al) = B/cpo, 

where n is the average density of bosons (and FLs) and A is the area of the sample. We can 
write the London free-energy functional in a dimensionless form as follows: 



(2.5) 



where 



A = ^ = ^^, ;3 = /l = H££^, A = A. (2.6) 

aQQ^m aQ^/2eieQ kTb kT oq 

All lengths are measured in units of ao and energies in units of for bosons and eo^o for 
FLs. We next discretize the integral along the 2;-axis by dividing it into M segments: 



■ ^ / p ^ m=l 1=1 



M N 

-3[R,,„]/fcT 



(2.7) 



where r = 13/M, m labels the planes and 



kT 



E 



■i,m+l 



R, 



2Ah 



(2.8) 



i<j 



We work only in the limit where f3 is large, which amounts to taking Lz very large and in 
the mapping onto 2D bosons corresponds to the ground state of the bosons at the absolute 
zero temperature (Tb = 0). We used /? = 375 and discretized the 2;-axis into M = 75 planes. 

We used the matrix-squaring method^^^ to calculate the right action so that we can 
work with a small number of planes along the 2;-direction. Working with the primitive action 
requires the use of a large number of slicing of the z-direction and is very time consuming^. 
The boundary conditions in the z-direction for a system of bosons are Ri(/3) = Rj(0), 
with all permutations j = P{i) of the indices being allowed. This is what is meant by 
the summation over P in Eq. ()2.7j) . For small g (a small repulsive interaction), which 
corresponds to large A, we expect a Bose-Einstein condensed phase where permutations are 
important. This is the superfluid phase. For small A, repulsion is large and permutations are 
rare as the Bose system is in its classical phase which is a Wigner crystal. These two phases 
correspond respectively to the FL liquid and solid phases. The melting line represented by 
A = Am is given approximately by the expression Bm = const. /T^. Near Tc this behavior is 
more complicated (see below). 

We now discuss the method of simulations of BSCCO : 

Because of the high anisotropy of the BSCCO system one can not use the simple picture 
given above. Here, instead, we follow a different model which can be cast in a form analogous 
to YBCO. We take the Lawrence-Doniach model^'' as our starting point for BSCCO. It leads 
to the following form of the London free-energy for inter-layer (IL) Josephson couplingi^^: 



,„,a„.)^i^(i..„(^ 



for |Ri_m - Ri,m+i| < 2r„ , and 



S/L(Ri 



^' 9 



|Ri,m Ri,m+l| 



(2.9) 



(2.10) 



for |Rj,m — Ri,m+i| > , where d is the inter-layer spacing and Vg is the healing length 
defined by Vg = '-yd. For the in-plane (IP) coupling we use: 



S7p(Rjj. 



I ^ij.m 



(2.11) 
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In principle one should add to the above interaction an electromagnetic interaction among 
the pancake vortices^-. The electromagnetic interaction becomes dominant in the limit of 



infinite anisotropy (7 00). Ryu et a/.— argue, (see also Ref. ISOl ) that for the value 
7 = 50 for BSCCO the Josephson coupling still dominates by one order of magnitude over 
the electromagnetic coupling. Their argument goes as follows: Clem''^ shows that if one has 
a straight array of pancake vortices along the 2;-axis, and one pancake is displaced a distance 
R in the x-direction than the magnetic energy of the configuration increases by an amount 



A^(«) = S^(C + .n^-J+A'„^-JJ. (2.12) 

where C is Euler's constant (=0.5772...). For large R {R^ A), the modified Bessel function 
Kq decays exponentially and thus the energy increases like ln(i?/A). For small R the Bessel 
function can be expanded in a power series in R/ A 

Ko{R/X) = - ln(i?/2A)(l + R^AX^ + ■■■)-€ + R\1 - C)/4X^ + ■■■ , (2.13) 

and thus the magnetic energy behaves like R^ to leading order in R which is the same as the 
quadratic behavior of the Josephson energy in Eq. ()2.9|) above. The ratio of the coefficients 
of the quadratic terms in the magnetic and Josephson energies goes roughly like 7^((i/A)^ 
(where (i/A ~ 1/100 for BSCCO). Thus for anisotropy 7 = 50 we get a factor of 0.25 (a 
somewhat more precise estimate^^ gives a ratio of about 0.1). Thus the magnetic interaction 
is negligible compared to the Josephson interaction for 7 = 50. For samples with 7 = 200 
these interactions are already comparable. For large values of R the magnetic interaction 
increases logarithmically and the Josephson interaction increases linearly so the magnetic 
interaction is always negligible. The key to the estimate given above is to consider not just 
two pancake vortices but a whole line with one displaced pancake. This argument is valid 
if the deviations of the vortices from straight lines are not too large. As for the in-plane 
interaction Clem showed that a linear array of pancake vortices gives rise to exactly the 
same magnetic field at a distance R away from it as produced by an Abrikosov vortex line. 
Thus Eq. fl2.11|) is consistent with the magnetic interaction of pancake vortices, again when 
the FLs do not deviate too much from straight lines. 

Equations (j2.9|) - (j2.1ip for BSCCO can be cast in a form similar to that for the YBCO 
with the following substitutions: 

5oao ci V 2(1 + M^))' kT- ^ ' ' 



With these changes the London free-energy functional would look like 



for (|Ri,^ 



kT 

R.i,m+i|) < 2rg, and 



2AV 



(2.15) 



S/L(R'i,r 
kT 



fR,; 



R 



i,m+l ) 



2Ah 



(2r, - (|R, 



R 



'i,m+l 



1))^ 



2Xh 



(2.16) 



for (|Ri,m — Ri,m+i|) > 2rg, when now again all lengths are measured in units of oq. While 
doing simulations at a fixed magnetic field B and temperature T , the term r|/2A^r will 
remain constant and would drop out of Ai? term in the Boltzmann factor. It however needs 
be considered during the measurement of the energy. The second term in Eq. ()2.1(jj) can be 
easily handled at the last stage of the Bisection method (see the next section for a discussion 
of the Bisection method). 

We can make use of a reduced temperature variable to make some expressions look 
simpler. First, using the fact that the temperature dependence of A arises mainly through 
Eq and neglecting the logarithmic corrections, one geta^^^^ 

Ao 



A 



1 



£0 OC — OC 



T 



and hence 



A oc 



Defining reduced temperature as 



T 



1 - T/T, 



we obtain 

A oc Tr\fB. 

This shows that the equation for the melting line is approximately 

5™ = const.(l-T/T,)V^'- 



(2.17) 



(2.18) 



(2.19) 



(2.20) 



Note that some authors^ use a temperature dependence of 1/a/1 — T"^ /T"^ in A, or eveit^ 



1/a/1 — T'^/T^. All these choices coincide near T^. The choice of temperature dependence 
of A is not expected to have a significant effect on the results. 
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III. NOTES ON THE SIMULATIONS 



The technique that we use to simulate our system is called Multilevel Monte Carlo Sim- 
ulation (MMC)^. There are several advantages in using this technique for the simulation 
of the FLs over the usual metropolis Monte Carlo (MC) method. In the discrete model, we 
work with FLs with the 2;-axis discretized into M planes, thus resulting in beads in 
each of the M planes. In the usual MC method one would displace a few of these beads in a 
plane by small random displacements inside a two dimensional box and then would accept 
or reject the move based on a probability given by the Boltzmann factor. 

A big disadvantage of using this technique is that it is difficult to move beads appreciably 
from their original positions over a number of MC steps. The reason for this is that a bead in 
a plane belonging to a PL finds itself in a local harmonic potential generated by the kinetic 
energy term involving this bead and the beads belonging to the same FL on either side of 
the plane. This harmonic potential becomes stronger and stronger at lower temperatures 
and magnetic fields. As a result, in the usual MC simulations beads keep moving around 
inside these local harmonic cages and end up sampling only a small part of the phase 
space. The other problem with the usual MC method is that there is no natural easy way 
of implementing FL cutting. If there are two FLs twisted around each other and if it is 
energetically favorable for them to reconnect each other in such a way as to lower their 
free energy then this step should be permitted in the MC method without regard to the 
question if this process occurs in reality. This is so because in the MC simulations phase 
space is sampled according to the probability distribution and all one needs is to generate 
configurations weighted by the Boltzmann factor, and the path followed in configuration 
space has nothing to do with any real dynamics. 

These two main drawbacks of the usual MC method are easily overcome in the MMC 
technique. First, one moves bigger chunks of FLs encompassing beads in several planes. This 
way one can avoid local harmonic traps. (This is like taking an aerial route to a destination 
rather than going through the zigzag maze of roads.) This is much in the spirit of Fourier 
space Monte Carlo where one first samples modes with smaller wave numbers and then move 
onto higher modes. 

The method of creating new FL configurations is based on the concept of the conditional 
probabilities. It is called the Bisection method^ because one starts sampling beads by 
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iteratively bisecting the Fls. At each stage of the division, the beads belonging to that 
stage are moved with some conditional probability factor Pj. It is important to make sure 
that the probabilities are chosen in such a way that detailed balance is satisfied at each 
stage of the division. One notes that the Pj's may not be the actual Boltzmann factors for 
the beads to be moved at different levels. But what is required is that when all Pj's are 
multiplied together, they cancel in such a way so as to leave the correct Boltzmann factor 
for the whole move. Thus, inherent in this algorithm is the fact that a move would finally 
be accepted only if it has been accepted at each stage of the Bisection method. The power 
of this method lies in choosing the appropriate Pj's. If these conditional probability factors 
are chosen judiciously, most of the rejections would take place at the initial stages of the 
bisection process when not too much computational effort has been spent yet. 

The cutting and reconnection of FLs is implemented naturally in a MMC method: per- 
mutation among the 3 or 4 lines chosen to be moved becomes the first among the many 
hierarchical steps one goes through before a move is finally accepted and the position of the 
beads updated accordingly. We typically moved a total of 15 to 20 beads distributed over 
5 planes. Permutations were sampled by a random walk algorithm through the space of 
permutations^^ (see Appendix B). 

In the case of YBCO we worked with a field of P = 4000 G. Working in the primitive 
approximation of the action would require the use of a smaller value of the dimensionless 
parameter r, which would require slicing of the z-direction into a large number of planes. 
To avoid this, the matrix-squaring method^^^ has been used to get the effective action for 
bigger values of r. For example Nordborg and Blattep^ use a value of r = 3.0 and they 
work with 100 planes. In this simulation a value of r = 5.0 has been used and the 2;-axis 
has been sliced into 75 planes. Choosing a bigger value of r by utilizing the matrix-squaring 
method makes it easier to equilibrate the system as compared to working with the primitive 
approximation. 

For BSCCO, we did not use the matrix-squaring method because of the complications 
involved due to the few extra terms in free energy which contribute depending on whether 
Rij is smaller or bigger than r^. Here, we used the natural inter-layer spacing d to calculate 
the parameter r at different temperatures and then used the MMC technique to efficiently 
sample the configurational space. 

As mentioned previously, in the present simulations we included only the Josephson 
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coupling. This approximation works well with YBCO but it could be less satisfactory for 
BSCCO because of its high anisotropy. For very anisotropic materials the electromagnetic 
coupling becomes import an t As discussed in Section 2, Ryu et al}- estimated that for 
anisotropy of magnitude 7 = 50 the Josephson interaction still dominates over the electro- 
magnetic interaction by an order of magnitude, but this will not be the case for 7 = 200 
which can characterize some samples. Olson et alM discuss how to include the electromag- 
netic interaction in a MC simulation, but they only consider the opposite limit where the 
Josephson coupling is totally neglected. To our knowledge there is no satisfactory treatment 
of both couplings included on equal footings. We carried out preliminary simulations which 
show that the inclusion of the electromagnetic coupling does not shift the position of the 
transition line much at a field of 125 G, thus supporting our current conclusions. These 
results will be reported elsewhere^. 

Simulations were usually carried out for 36 and 64 FLs. The decay of the structure 
factors at the transition temperature becomes sharper when one uses a larger number of 
FLs. However, no appreciable shift in the transition temperature is seen while working with 
the smaller (A^ = 36) or the larger system (A^ = 64). We did not run our simulations for 
even larger systems as it becomes computationally very time consuming. Typical simulation 
times were ~ 3.5 x 10^ MC steps. Each MC step involved moving 3-4 lines in 5 planes 
simultaneously. We usually averaged over 10-15 realizations of the columnar disorders, 
though some results have an average over as many as 20 realizations of the disorder. 

Columnar disorder is modeled as an array of straight cylindrical wells of typical radius 
rrf=25-35 Aplaced randomly throughout the cross-sectional region of the sample and oriented 
along the 2;-directioni^. Each columnar defect is of length L^. The density of the columnar 
pins can be varied by changing their number for a given cross sectional area, and the strength 
is controlled by a positive dimensionless parameter rj. If a bead happens to wander inside 
a columnar well, we include an extra free energy of — 77^^^-. The defect concentration was 
taken to be a 20 percent ratio of defects to FLs which means B^/B = 0.2. The strength of 
disorder was set at 77 < 0.5 for BSCCO. For YBCO, 77 = 0.5 was found to be too large in the 
sense that the transition became too broad and useful information could not be extracted. 
Thus, for YBCO we kept r] < 0.3. 

Other parameters used for YBCO are 7 = 5, Aq = 1500 A, ^^^(O) = 15 A, = 12 A. 
Parameters for BSCCO are as follows: 7 = 50, Aq = 1700 A, ^„f,(0) = 20 A, = 15 A. 
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IV. MEASURED QUANTITIES 

In this section we describe many different physical quantities that were monitored during 
the simulation. From the variation of these quantities with the temperature we can extract 
important conclusions about the different phases of the system. 



A. Energy 



In terms of the reduced temperature, the energy 



can be simply written as^: 



E = T'-^HE{A,(3,N)) 



E = TriSi + S2] 



where, 



and 



2 _ ^2 
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2A2r 



for (|Ri,„ - Ri,m+i|) < 2r, 



5' 



(Ri.m Ri,m+l) 



2A2r 

for (|Ri,m - Ri,m+i|) > 2rg, whereas 



(2rg — (|Rj,m — R-i,m+l|))^ 



2A2 



r 



^2 = ( E ( 

m,i>j ^ 



A 



(4.1) 



(4.2) 



(4.3) 



(4.4) 



(4.5) 



Any discontinuous jump in E would indicate a FOT. From this discontinuous jump in energy, 
AE, we can also calculate the jump in entropy As, 

AE 



As 



T 



(4.6) 



B. Translational structure factor 

The translational structure factor S^Qi-^^i^) is defined as, 

SiQi^i.) = ^(EeOQ'-'-(^--^-))) (4.7) 

ij,m 
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where (...) stands for the MC average, Q/1,/2 is a reciprocal lattice vector and is of the 
form 

Qi.^i, = hh^ + /2b2 (4.8) 

and li and I2 are some integers, bi 2 represent the basis vectors of the reciprocal lattice 

2n 

bi,2 = r^(ei,2 - e2,icos6') (4.9) 

ao sm & 

where ^ = 7r/3, ao is the nearest neighbor distance and 61^2 are the unit vectors along the 
hexagonal unit cell such that 

ei ■ e2 = cos 6. (4.10) 

In the simulations, li and I2 in Eq. ()4.7|) are chosen to be 1 and or and 1. These 
choices correspond to the first Bragg peak. We will normally write 5'(Qi^o) as simply S'(Qi). 
This quantity gives important information about the phase of the system. In an ordered 
phase where Fls sit on a triangular lattice, >S'(Qi) is of the order of A^. In the disordered 
phase, it saturates to almost zero as the system size increases. In the simulations, however, 
there is a problem with measuring ^(Qi), especially in the presence of disorder. We will 
return to this point at the end of the section. 

C. Hexatic Structure Factor 

We use Delaunay triangulation to measure the hexatic order parameter, ip^, which is 
defined as, 

^6 = (EE— E^^^''^^^'"^)' (4.11) 



m=l 1=1 ' j=l 

where Zi^m denotes the number of the nearest neighbors of a bead at position i, m and it is 6 
for a perfect hexagonal lattice, Oij^m stands for the bond angle, that is the angle that vector 
Hij^m makes with an arbitrary axis. Just like the ^(Qi), this quantity too has a large value 
in an ordered phase and saturates at a finite value for a system of finite size. 

D. Line Entanglement 

As we allow permutations of FLs, we can define a number N^/N as that fraction of the 
total number of FLs which belong to loops that are bigger than the size of a "simple" loop. 
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A simple loop is defined as a set of M beads connected end to end, M being the total 
number of planes. Loops of size 2M, 3M... start proliferating at and above the transition 
temperature and in the corresponding 2D boson system this proliferation is related to the 
onset of the superfluidity. 

E. Line Wandering 

The transverse FL fluctuations are measured by 

u\z) = ((R(^o + z)- R(zo))')/a^, (4.12) 

which is independent of zq. At the transition temperature v?{z) undergoes a large increase 
for large z. 

We want to emphasize one important aspect of measuring the translational structure 
factor. The usual way of measuring this quantity is by choosing Q^^^^a corresponding to 
the first Bragg peak i.e., Qi in Eq. ()4.7|) . Now, it might happen that the configuration 
of the FLs comes close to making an almost perfect hexagonal lattice but its basis vectors 
are not aligned along the usual major axes of the rhombically shaped cell encompassing the 
system. If we use the reciprocal lattice vector Qi to measure the structure factor of such 
a configuration, we would end up getting a very small value for S'(Qi) and might wrongly 
conclude that the system is in a very disordered state. This happened many times in our 
simulations; we got a very low value of the translational structure factor while the hexatic 
order was indicating a high degree of orderliness in the FL lattice. To remedy this situation 
the translational structure factor is measured at 60 different angles, 1 degree apart, and 
choose that number which gives the largest possible value of the S'(Qi(a;)) corresponding to 
some angle a. After implementing this technique we find that >S'(Qi) and ip^, which differ 
so much initially, almost follow each other. 

V. RESULTS 

1. For YBCO, simulations were carried out at a magnetic field of -B = 4000 G. At this 
field we have 2ao ~ Aq. By the argument given in Ref. 5, we conclude that our results 
should hold qualitatively for any B such that < A. We checked our simulation results 
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FIG. 1: YBCO: Structure Factor at the first Bragg peak as a function of A at various disorder 
strengths (r/'s) for = 36. A clear shift in transition point towards higher A's is seen for r] > 0.2 

against those by Nordborg and Blatter^^. First we verified that our code was working fine 
by comparing our results against those given in Ref. Q. Working in the hmit of — oo, 
we kept r fixed as A was increased. We got a sharp transition at A = 0.0605. A jump was 
also seen in the energy as defined in the previous section and was found to be 0.0 13A;, again 
the same as in Ref. llj 

The effect of introducing columnar disorder on the structure factor at the first Bragg 
peak is shown in Fig. ^ In Fig. |21we depict the jump in the energy for a pure as well as for 
a disordered system. 

We note that at lower strengths of the disorder, the melting line is almost unaffected 
and the transition takes place at A 0.0605. However, as we increase the strength of 
the disorder, the melting transition shifts towards higher values of A which means that the 
melting line shifts towards higher temperatures and/or higher magnetic fields. For t] up 
to 0.1 no change is seen in the melting curve or in the jump in the free energy functional. 
However, at rj equal to 0.2 the structure factor comes down at around 0.0625 and, the jump 
in the energy becomes gradual. This finding is in agreement with several experimental 
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FIG. 2: YBCO: Energy as a function of A for A'^ = 36. A jump in the energy is seen up to r/ = 0.1. 
For r] > 0.2, the rise in the energy is soothed out. 

studies where a change in the irreversibihty hne is seen with the introduction of columnar 
disordei.2L2mim. 

Next, we look at the PL entanglement. A^^e denotes the number of FLs which do not form a 
simple loop. By looking at the N^/N vs. A graph in Fig. IHl it is clear that the entanglement 
is suppressed by as much as almost one order of magnitude just after and in the vicinity of 
the original transition line at A = 0.0605. This result is expected, as it is well known that 
point disorder helps in line wandering and entanglement while the columnar disorder has 
just the opposite effecti^^ since a nearby FL is induced to align along the columnar defect, 
and thus its transverse fluctuations are reduced. Unfortunately, there have been few studies 
on columnar disorder in YBCO. Sen et alM work with vanishingly small magnetic field and 
infinite pinning strengths, in the vicinity of the lower branch of the melting line. 

Figure |3] shows the effect of columnar disorder on a system of 64 FLs in YBCO. This 
graph shows a transition at almost the same A as in Fig. ^ This is a confirmation that 
finite size effects are not important in our simulations. 

2. For the BSCCO system the E vs. T graphs for different vj values are depicted in Fig. 
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FIG. 3: YBCO: Fraction of entangled FLs as a function of A. For rj < 0.1 a sharp jump in the FL 
entanglement is seen. For r] > 0.2, the rise of N^/N with A becomes gradual as compared with the 
clean system. 

0(36 FLs) and Fig. H (64 FLs). 

For the pure system (77 = 0) we see a sharp jump in the energy at exactly the point 
where the translational as well as hexatic structure factors sharply decline. This jump in 
energy is a clear signature of a first order transition. From the energy jump we can compute 
the jump in entropy (As). The jump in energy at the transition is around AE = 12.0k K 
per vortex per layer, which gives As = AE/Tm = 0.19fc which is small compared to the 
experimental value of As = 0.40A; at 5 ~ 125 G. However, these values for As, Tm and B 
are in qualitative agreement with the same system studied using a different mode^^. 

Introduction of the columnar defects of a finite strength shows some interesting effects. 
We put columnar pins at random positions with a concentration fixed at 20 percent of the 
FLs. We study the BSCCO system for pins of low strength [i] = 0.2) as well as for a higher 
strength (77 = 0.5). Columnar disorder of strength up to 77 = 0.2 appear to have little effect 
on the system . This can be seen from the energy vs. temperature graph in Figs. El and El 
where the curve for the pure system and the one in the presence of low disorder strength 
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FIG. 4: YBCO: Structure Factor at the first Bragg peak vs. A for = 64. Smootli lines are 
provided as a visual aid. Transition for rj = 0.0 and ij = 0.2 are seen almost at the same A's as 
with N = 36. 

[r] = 0.2) fall on top of each other, except in the transition region where the jump in case of 
the system with columnar pins is more gradual. Also, from Figs. [7| and |H1 we see that the 
S'(Qi) is not much affected in the presence of the columnar disorder of strength rj = 0.2. The 
same is true for the hexatic order depicted in Fig. El This means that for the lower strengths 
of the disorder, the vortex-vortex interaction is dominant in the range of the temperature 
where we carried out simulations. As a result the first order vortex-lattice melting transition 
(FOT) line is almost not affected. 

The nature of FL melting changes, however, when the strength of the disorder is in- 
creased to T] = 0.5 as seen in Figs. [7| and |H1 We see that at the lower temperatures the 
order parameter is suppressed. As we increase the temperature, the order parameter starts 
to rise and joins the melting curve of the pure system and then falls along it at even higher 
temperatures. Fig. El indicates that the hexagonal structure factor also shows a similar be- 
havior. This clearly shows that the FLs start to disengage from the pinning centers as they 
wiggle more. At higher temperatures, columnar pins with the strength chosen here are not 
effective in pinning the vortex system. This becomes clear from the fact that the melting 
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FIG. 5: BSCCO: Energy as a function of the reduced temperature for rj = 0,0.2,0.5 at B = 125 
G for = 36. Lines were added to help visuahze the energy jumps. A discontinuous change in 
energy is seen at Tr = 250, 260 K (T ~ 66.2, 66.9 K) for rj = 0, 0.2 respectively. For r] = 0.5 only a 
change in slope can be seen at Tr ~ 268K (T ~ 67.4 K) instead of a discontinuous jump. 

curve for the pure system and the one with columnar disorder join each other at a tempera- 
ture below the FOT. Further interpretation and a possible explanation of this phenomenon 
will be discussed in more detail in the Conclusions section below. This convergence of the 
curves for the pure and disordered cases is also borne out in the E vs. T graphs in Fig. El 
and Fig. IHl Initially at low temperatures there is a big difference in energies of the systems 
with no disorder {r] = 0.0) and high disorder [t] = 0.5). However as the temperature is 
increased, the two curves come closer and finally merge together in the liquid phase. The 
jump in energy in presence of columnar defects can still be seen in the E vs. T graph for 
a disorder strength of ?7 = 0.2. This tends to suggest that the FOT in a BSCCO system is 
not affected by disorders kept at as many as 20 percent of the total number of lines used in 
the simulations as long as we keep the strength of the disorder less or equal to rj = 0.2. For 
rj = 0.5, the rise of energy with temperature becomes much more gradual and we do not see 
any discontinuous jump in energy at any temperature. On the other hand an abrupt change 
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FIG. 6: BSCCO: Energy as a function of the reduced temperature for rj = 0,0.2,0.5 at B = 125 
G ioi N = 64. A jump in energy can be seen at « 238, 242 K (T w 65.3, 65.6 K) for r? = 0, 0.2 
respectively. Again no discontinuity is observed for r/ = 0.5, but only a change in slope at Tr ~ 255 
K (T « 66.5 K). 

in slope of the energy vs temperature graph is observed, suggesting a discontinuity in the 
specific heat characterizing a second order transition. 

In Fig. Uniand Fig. [TTl the snapshots of FLs, projected on a plane, at temperatures less 
than the transition temperature and at a temperature bigger than the transition temperature 
are shown. From Fig. it is easily seen that at low temperatures most columnar defects 
have captured FLs. Also, the FLs make simple loops and have cleverly set themselves so 
as to make a hexagonal lattice and yet occupy as many defects as possible. The transverse 
fluctuations of the trapped FLs are greatly suppressed at low temperatures. At higher 
temperature beyond the transition point, we see that columnar defects are not occupied any 
more and a lot of FLs are entangled. 

Inspection of snapshots like Fig. [TUl gives support to the assertion of Sen et al. that 
the Bose glass consists of patches of ordered regions with only short range positional and 
orientational order. This phase is different from the Bragg-glass in systems with point pins 
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FIG. 7: BSCCO: Structure factor at the first Bragg peak vs. temperature at i? = 125 G for 
= 36. No shift in transition point is seen for 77 = 0.2. In the presence of columnar pins with 
f] = 0.5, the transition temperature increases from 66 K to 68 K. Also, we can see that 5'(Qi) 
starts to rise close to 60 K 

which is characterized by quasi-long-range order. 

Figure IT^ shows the mean squared displacements of the FLs in the z-direction at different 
temperatures for a pure system. At lower temperatures u^{z) saturates for large z but in 
the liquid state it grows linearly^^. A large gap in v?'{z) for large values of z seen at the 
transition temperature signals the onset of the entanglement of the FLs. 

In the presence of the columnar disorder of the strength 77 = 0.5 we see from Fig. 
that u^{z) at temperatures less than the transition temperature is suppressed compared to 
the corresponding u'^{z) in the pure system. This result is in agreement with the findings in 
Ref. [3 Also the big gap in v?[z) occurring at the melting transition has moved towards a 
higher temperature signaling a shift in the position of the melting transition in the presence 
of columnar defects. 
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FIG. 8: BSCCO: Structure factor at the first Bragg peak vs. temperature at B = 125 G for 
N = 64. For r] = 0.2 no significant change in the structure factor is observed. In the presence of 
columnar pins with rj = 0.5, transition temperature increasing from 65 K to 67 K. Also, we can 
see that 5(Qi) starts to rise close to 61 K 

VI. CONCLUSIONS 

For lower strengths of disorder [i] < 0.1 for YBCO and rj < 0.2 for BSCCO) no appreciable 
shift in the melting line was seen. For these strengths, a sharp drop in the translational 
and hexatic structure factors takes place at the transition. Also, the jump in energy at 
the transition is not affected much as compared with the case when there is no disorder 
present. This suggests that for smaller concentrations of the columnar defects the transition 
still remains first order. This result is in agreement with Ref. where even with the 
introduction of the columnar disorder no shift in transition temperatures is seen as long as 
the concentration of columnar disorder introduced is small. 

It is found that for YBCO as well as BSCCO, the melting transition shifts towards higher 
values of temperature and magnetic field when random disorder is introduced provided its 
strength exceeds a threshold which is different for YBCO and BSCCO. The size of the shift, 
for a given concentration of defects, depends on the strength of the disorder. For YBCO, 
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FIG. 9: BSCCO: Hexatic order vs. temperature at B = 125 G for iV = 64. When compared with 
the previous figure, it can be seen that ipe almost follows S'(Qi). 
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FIG. 10: BSCCO: A typical configuration in the solid phase (low temperature) for A'' = 64 FLs for 
B(j)/B = 0.2. FLs have been projected onto one plane. 
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FIG. 11: BSCCO: A typical configuration of FLs in the liquid phase (high temperature) for A'^ = 64 
and B(f,/B = 0.2. FLs have been projected onto a single plane. Columnar defects are not drawn to 
scale. Some FLs on the boundary do not seem to make loops. That is only because virtual images 
of FLs outside the cell are not shown. 

a considerable shift in the melting line towards higher temperatures and magnetic fields 
was seen for r] = 0.2 and rj = 0.3. This shift was bigger for rj = 0.3 than for rj = 0.2. 
Similarly, for BSCCO, a large shift in the melting line towards a higher temperature was 
found for t] = 0.5. These results are in tune with numerous experimental finding a^^i^^i??!?? as 
well as theoretical predictioit^^ where the irreversibility line is seen to shift towards higher 
temperature and magnetic fields in the presence of columnar defects. The qualitative reason 
for this effect is that due to the interaction with the columnar defects the transverse thermal 
fluctuations of the FLs are reduced and thus the melting transition, as determined from the 
Lindemann criterion, takes place at at a higher temperature. 

At 1] = 0.5 the jump in energy is not discernible any more. Instead, a change of slope 
corresponding to a specific heat discontinuity is observed. This means that the transition is 
probably not of first order in the presence of columnar disorder of higher strengths, but it 
is rather a continuous (second order) transition. 
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FIG. 12: BSCCO: Line wandering along the z-direction at B = 125 G for AT = 64. Here d is the 
distance between adjacent planes. A large increase in line wandering occurs at T 65 K. 




FIG. 13: BSCCO: Line wandering along the z-dircction at B = 125 G for = 64 in presence of 
disorder of strength rj = 0.5. The big jump in u'^{z) has moved to T = 66 K now. 
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The most dramatic outcome of this study for BSCCO is that for some values of the apphed 
field and defect strength both the translational and hexatic structure factors start to rise at 
a certain temperature as the transition is approached from the lower temperature side. This 
is an unusual result, that to our knowledge, has not been seen in previous simulations. This 
fact can be explained as follows. At low temperatures, the free energy is dominated by energy 
effects rather than entropy considerations, and pinning effects are dominant. As a result 
most columnar defects capture a FL, while the rest of the FLs adjust themselves in positions 
such as to minimize the free energy. However, as the temperature is increased, FLs start 
to decouple from the defect sites, which allow the vortices pinned at interstitial positions 
to move themselves into a more ordered arrangement, thus increasing the structure factors 
compared to the situation at lower temperatures. As was mentioned before, due to its high 
anisotropy, FLs in BSCCO behave more like a collection of two dimensional pancakes rather 
than rigid rods. These pancakes are comparatively more difficult to get pinned all at once 
by a columnar defect. On the other hand, FLs are much stiffer in YBCO and it is easier for 
a columnar defect to capture a FL all along its length. It can be shown (Eq. (9.49) in Ref. 
Q) that the depinning temperature for FLs, Tdp, is inversely proportional to the anisotropy 
parameter 7. Thus, we expect FL depinning to occur at comparatively smaller temperatures 
in BSCCO than in YBCO. Since the melting temperatures are not that different between 
these materials at the fields we consider, the depinning process for BSCCO occurs further 
below the melting temperature than in YBCO. This allows the structure factor to increase 
above the depinning temperature before its ultimate decline at the melting transition. For 
YBCO, the depinning takes place close to the melting temperature, and it is difficult to 
detect any rise in the structure factors, especially for small systems, because it is masked by 
the decrease in order due to the stronger thermal fluctuations. 

In principle. The rise of the order parameters in the presence of low amount of columnar 
defects as the melting transition temperature is approached could be observed experimen- 
tally, if the appropriate parameters are tuned correctly. In small angle neutron scattering 
(SANS)^^^ one can measure the integrated intensity over a Bragg peak of wave vector Qi 
which is proportional to the translational structure factor ^(Qi) measured in the simulation. 
Another commonly used technique is muon spin rotation^^ which gives information about 
the width of the magnetic field distribution in the sample. This is not directly proportional 
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to the structure factor at a given Qi, but is rather given by 




(6.1) 



where {u^) is the mean square deviation of vortices from their average position. It is possible 
to measure this quantity in the simulations but this was not done it in the present work. 
This quantity might not show the unusual rise describe above since it is dominated by 
(v^) which is very likely monotonically increasing with temperature yilding a monotonically 
decreasing line width. Thus in order to look for the effect observed in this paper we suggest 
using the SANS technique to look at a BSCCO sample with columnar pins described by 
a matching field of about 25 G and an applied field of about 125 G. It is not clear to us 
what is the corresponding r] parameter describing the pinning strength of the experimental 
defects. According to Blatter et al^ rj lies in the range 0.1-1. 
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APPENDIX A: ENERGY SUM OVER THE IMAGES 

Consider a rhombically shaped region with side L and angle 9, unit vectors are ei, e2, 
with ei • e2 = cosO. In practice we took 9 = 60° but we leave the discussion here more 
general. The Green's function Go which describe the 2D coulomb interaction between one 
vortex and another including all its images, as is implied by the periodic boundary conditions 
is given by the solution to London's equation (see e.g. Ref. l| 



with the parameter A setting the scale for the range of the interaction. The solution is given 



(1-A2v')Go(R,A) = A25(R), 



(Al) 



by 




(A2) 



with 



R = i?iei + i?2e2, Q = riihi + n2b2, 



(A3) 
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where Q runs over all reciprocal lattice vectors spanned by 

27r 



Bj — ej cos 6), (A4) 



' Lsin^^ 

for (ij) = (1,2), (2, 1) . Substituting in Eq. (jMj) we obtain 

r rR AW ^ V V exp[t^{niRi + n2R2)] 

ni=— oo ri2=— oo / \ / i z 

We are now going to carry out the summation over rii analytically. This can be done by 
using the formula 

oo 

f{n) = — (residues of 7r/(z)(cot ttz — i) at poles of f{z)). (A6) 

n=— oo 

We have subtracted the constant i from cot ttz to ensure that \zf{z) (cot ttz — ^ on the 
contour of integration when z has a large negative imaginary component. The contour of 
integration is a square with sides parallel to the real and imaginary axes with the origin in 
the middle, in the limit that its size goes to infinity. In our case 

^(^) = (z-/3-z7)(z-/3 + z7)' ^^^^ 

with appropriate values of /? , 7 , and x . This function has two simple poles, and the 
residues can be easily evaluated. The final answer becomes (relabeling as ) 

r (-R \^ - ^^^^ cos{t2n - 2Tipn) sinh(7„ti) + cos(t2n) sinh(7„(27r - ti)) 

Uo[tt,A)- 2 ^2_.^ 7„(cosh(2vr70 - cos(2vr/?0) ' ^"^^^ 

where 

27ri?i 2tt 



ti = —^, t2 = —{Ricose + R2), (3n = ncose, 7„ = sin^Vn2 + LV(27rA)2. (A9) 

This expression is simpler than the one used by Nordborg and Blatter since it does not have 
different expressions for even and odd n. It also converges faster in certain regions. 



APPENDIX B: PERMUTATION SAMPLING 

Essentially the same method is used for permutation sampling as was used in Ref. Iiol 
The only difference is that we use permutation space of only the neighboring lines. This 
is so because even if we were able to get a permutation step, involving a large number of 
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lines, accepted at the first stage of tlie algoritlim, it would be very likely to get rejected at 
the following stages. So we work with only 3-5 lines. Sufficiently long segments (typically 5 
planes) of a number of lines were cut. Care should be taken to make sure that chosen FLs 
are the nearest neighbors in the plane where the reconnection of FLs is going to take place. 
These points can be implemented easily with the concept of linked lists and pointers^. 
Also, care is to be taken that even though a PL may be far away from some other FL in 
the rhombically shaped unit cell, it can still permute with it through one of the images of 
the latter. These few simple points are very important to implement the whole procedure 
correctly. Just for a check, we tried with the sampling procedure given in Ref. H This 
gave results in good agreement with the sampling procedure given above. 
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